R ist eine öffentliche Statistik-Programmiersprache, welche sich für eine breite Anwendung einsetzen lässt.
RStudio ist in R-Kreisen die bekannteste Programmier-Entwicklungsoberfläche. Sie wurde von der Firma RSTudio entwickelt und es gibt eine open-saurce Version.
Die katholischen Kirche möchte wissen wie gross der Anteil der ständigen Wohnbevökerung an über 80 Jährigen in den Kirchgemeinden ist.
Inputdaten: Kirchengemeinde (shp); Stand 2019 städnigen Wohnbervölkerung (sph); Stand 2020
Methodik
Kirchengemeinde: Stand 2019 => Stand 2020:
Die Kirchgemeinde “Hirzel-Schönenberg-Hütten” wurde getrennt:
=> Hirzel hat sich der Kirchgemeinde Horgen angeschlossen.
=> Schöneberg und Hütten haben sich der Kirchgemeinde Wädenswil angeschlossen.
Die Daten der Bevölkerungsstatistik werden auf die neuen Kirchgemeinden aggregiert.
Hinweis:
Da die Bevölkerungsstatistik auf eine Hektare aggregiert werden, stimmen die aggregierten Zahlen nicht ganz mit der realität überein. Wir verwenden die Bevölkerungsstatistik da die Daten OGD sind und ihr so das Beispiel durchspielen könnt. Wir haben die Daten eigentlich auf EGID genau und können auf dieser Basis auswertungen machen.
Es werden alle packages geladen, welche für die Aufbereitung verwendet werden. Zudem wird eine Funktion erstellt, welche die Shapefiles auf nicht valide Geometrien untersucht.
# Laden aller benötigten Packages
library(sf) # Vektor-Library
library(dplyr) # Library zur Datenaufbereitung
library(stringr) # Library für regex
library(ggplot2) # Library zur Plot-Erstellung
library(ggspatial) # Library mit Map-Features (Nordpfeil usw.)
library(mapview) # Library für interaktive Karten
# Hilfsfunktion um die Validität der Geometrien zu prüfen. Falls eine Geometrie kaputt ist, gibt es eine Warnung. Die Topologie kann jedoch nicht geprüft werden.
check_validity <- function(sf_df){
if(FALSE %in% sf::st_is_valid(sf_df)){
warning("There is an invalid geometry in the Dataset. Use st_make_valid() to repair the geometry.")
}else{
print("The geometry is valid.")
}
}Der Import der Datensätze, in unserem Fall shapefiles, kann mit st_read() gemacht werden. Es können aber auch zum Beispiel gpkg, geojson und Daten aus einer DB eingelesen werden. Daten können auch, mit hilfe eines SQL Queries, partiell eingelesen werden.
# Importieren des Kirchgemeinden-Shapefiles des Jahres 2019
kg <- sf::st_read(paste0("../shapefiles/kirchgemeinden_2019.shp"), stringsAsFactors = FALSE, crs=2056)## Reading layer `kirchgemeinden_2019' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/kirchgemeinden_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 62 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2669245 ymin: 1223895 xmax: 2716901 ymax: 1283343
## Projected CRS: CH1903+ / LV95
# Importieren des Gemeindegrenzen-Shapefiles des Jahres 2020
gem <- sf::st_read(paste0("../shapefiles/gemeinden_2020.shp"), stringsAsFactors = FALSE, crs=2056, query = "SELECT GDE_ID, ART_N FROM \"gemeinden_2020\" ")## Reading layer `gemeinden_2020' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/gemeinden_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 162 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2669245 ymin: 1223896 xmax: 2716900 ymax: 1283343
## Projected CRS: CH1903+ / LV95
Das sf package hat eine eigene plot-Funktion welche fähig ist Karten zu plotten. Diese ist sehr einfach, es hilft aber um sich einen schnellen Überblich zu verschaffen. Mit mapview() gibt es eine Funktion, die erlaubt, eine Karte interaktiv anzuschauen. Dazu wird im backend leaflet verwendet.
# statischer Plot
plot(kg$geometry)# dynamischer Plot
mapview(kg)Das sf-package bietet Funktionen an um die Geometrie jedes polygons zu prüfen und zu korrigieren.
check_validity(kg)## [1] "The geometry is valid."
check_validity(gem)## [1] "The geometry is valid."
Für das Splitting der Kirchgemeinde, werden die politischen Gemeindegrenzen verwendet. Also der Teil der Kirchgemeinde, welcher in der politischen Gemeinde Horgen liegt, wird der Kirchgemeinde Horgen zugewiesen, der andere Teil der Kirchgemeinde Wädenswil. Um den Prozess relativ simple zu halten, extrahiere ich nur die benötigten Polygone.
# Filtern der Kirchgemeinde mit der Nummer 19 (Hirzel-Schönenberg-Hütten)
kg_19 <- kg %>% filter(Kirchnumme == 19)
# Filtern der politischen Gemeinden (Horgen, Wädenswil)
gem_h_w <- gem %>% filter(GDE_ID %in% c(293, 295))# interaktive Darstellung
mapview(gem_h_w, zcol = "ART_N") + mapview(kg_19, zcol = "ART_TEXT", col.regions = c("red"))Um dann die gesplitteten Polygone den richtigen Kirchgemeinden zuordnen zu können, wird eine Mapping-Tabelle erstellt.
# Erstellung einer Tabelle mit tibble
kg_fus_tibble <- tibble(
GDE_ID = c(293, 295),
Kirchnumme = c(45, 21),
GEMEINDENA = c("Wädenswil", "Horgen"),
ART_TEXT = c("Kirchgemeinde", "Kirchgemeinde"),
ART_CODE = c(1,1)
)
# Darstellen der Tabelle
knitr::kable(kg_fus_tibble)| GDE_ID | Kirchnumme | GEMEINDENA | ART_TEXT | ART_CODE |
|---|---|---|---|---|
| 293 | 45 | Wädenswil | Kirchgemeinde | 1 |
| 295 | 21 | Horgen | Kirchgemeinde | 1 |
Mit dem Verschnitt der beiden Layers wird die Kirchgemeinde in mehrere Polygone unterteilt, welche die ID der darunterliegenden politischen Gemeinden enthalten.
# Verschneiden der Kirchgemeinde mit den politischen Gemeinde um die Kirchgemeinde zu trennen.
#Der Output wird als Multipolygon ausgegeben
union_gem <- st_intersection(kg_19, gem_h_w) %>% st_cast("MULTIPOLYGON")
mapview(union_gem, zcol = "GDE_ID")# Multipolygon -> single Polygon
union_single <- st_cast(union_gem, "POLYGON")
mapview(union_single, zcol = "GDE_ID")Es gibt einige Sliverpolygone, welche den zwei grossen Polygonen zugeordnet werden.
# Die sliverpolygon kleiner als 1000m2 werden der Kirchgemeinde Horgen zugeordnet.
# Wir können das so machen, da es nur sliver-polygone an der Grenze zwischen der
# Kirchgemeinde Horgen und der politischen Gemeinde Wädenswil auftreten.
union_merged <- union_single %>% mutate(area = as.numeric(st_area(.))) %>%
mutate(GDE_ID = ifelse(area < 1000, 295, GDE_ID)) %>%
group_by(GDE_ID) %>%
summarize()
mapview(union_merged, zcol = "GDE_ID")Die neuen Polygone enthalten zur Zeit nur die politische Gemeinde-ID und die alte Kirchennummer. Mit der Mapping-Tabelle werden nun die neuen Kirchgemeinde-ID hinzugefügt.
# Die mapping Tabelle wird nun an die beiden neuen Polygone hinzugefügt um die Kirchennummer den Polygone zuzuordnen.
new_kirchid <- union_merged %>%
left_join(kg_fus_tibble, by = "GDE_ID") %>%
select(-GDE_ID)Wir ersetzten nun die Kirchgemeinde mit der Nummer 19 durch die zwei neuen Polygone. Die neuen Polygone werden dann mit den bestehenden Kirchgemeinden gemerged um einheitliche Polygone zu erhalten.
# Ersetzen der alten Kirchgemeinde durch die zwei neuen Polygone.
# Dann werden die Kirchgemeinden zusammengefügt
kg_2020 <- kg %>%
filter(Kirchnumme != 19 | is.na(Kirchnumme)) %>%
bind_rows(new_kirchid) %>%
st_cast("MULTIPOLYGON") %>%
mutate(STICHTAG = "2021-01-01") %>%
group_by(ART_TEXT, ART_CODE, STICHTAG, GEMEINDENA, Kirchnumme) %>%
summarize()
check_validity(kg_2020)## [1] "The geometry is valid."
plot(kg$geometry)
plot(kg_2020$geometry)So, jetzt haben wir den Kirchgemeindestand vom Jahr 2020.
Um nun den Anteil der über 80 Jährigen zu berechnen, verwenden wir den OGD-Datensatz der Bevölkerungsstatistik. Da dies Hektardaten sind, ist es lediglich eine Annäherung. In unserer DB hätten wir Gbäude-Spezifische Daten, um das Beispiel nachspielen zu können, wird hier aber der OGD-Datensatz verwendet.
Da bereits die Anteile auf Hektarraster-Ebene berechnet wurden, müssen wir zuerst die Prozente wieder in Absolute Zahlen umrechnen.
# laden der Bevölkerungsstatistik
bevstatistik <- sf::st_read("../shapefiles/bevstatistik_punkte.shp")## Reading layer `bevstatistik_punkte' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/bevstatistik_punkte.shp' using driver `ESRI Shapefile'
## Simple feature collection with 36256 features and 23 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 2669650 ymin: 1224850 xmax: 2716150 ymax: 1283350
## Projected CRS: CH1903+ / LV95
# Vorbereiten der Daten
bevstatistik_abs <- bevstatistik %>%
mutate_at(vars(matches("_P")), ~round(PERS_N/100*.)) %>%
rename_at(vars(matches("_P")), ~str_replace_all(., "_P", "_A")) %>%
mutate_at(vars(-geometry), ~ifelse(. == 9980 | . == -999, NA, .))Die Punktdaten und Kirchgemeinden werden verschnitten um dann über die Kirchgemeinde aggregieren zu können.
# Verschneiden der Punktdaten mit den Kirchgemeinden
bev_within_kg <- sf::st_intersection(bevstatistik_abs, kg_2020)
mapview(bev_within_kg %>% filter(Kirchnumme==45))Die Punktdaten besitzen nun die Information der Kirchgemeinde und wir können die Daten nach dieser Variable aggregieren. Dann müssen wir noch den Anteil der über 80 Jährigen berechnen.
# aggregieren der Daten
anteil_ue80_pro_kg <- bev_within_kg %>%
st_drop_geometry() %>%
filter(!is.na(Kirchnumme)) %>%
group_by(Kirchnumme) %>%
summarize(anzahl_ue_80 = sum(J_80PLUS_A, na.rm = T),
anzahl_tot = sum(PERS_N, na.rm = T)) %>%
ungroup() %>%
mutate(anteil_ue_80 = round(anzahl_ue_80/anzahl_tot*100, 1))
knitr::kable(head(anteil_ue80_pro_kg, 10))| Kirchnumme | anzahl_ue_80 | anzahl_tot | anteil_ue_80 |
|---|---|---|---|
| 1 | 1136 | 19016 | 6.0 |
| 2 | 1143 | 26075 | 4.4 |
| 3 | 1148 | 22712 | 5.1 |
| 4 | 524 | 11770 | 4.5 |
| 5 | 795 | 13242 | 6.0 |
| 6 | 600 | 14330 | 4.2 |
| 7 | 1586 | 35178 | 4.5 |
| 8 | 1539 | 42131 | 3.7 |
| 9 | 1370 | 28268 | 4.8 |
| 10 | 2236 | 43619 | 5.1 |
Die Aggregierten Daten können nun dem Kirchgemeindenshapefile angefügt werden und wir können nun eine Choroplethenkarte erstellen.
kg_2020_with_anteil <- kg_2020 %>%
left_join(anteil_ue80_pro_kg, by = "Kirchnumme") %>%
select(Kirchnumme, anteil_ue_80)
# plot der Daten
ggplot(kg_2020_with_anteil) +
geom_sf(aes(fill = anteil_ue_80), color = "white") +
theme_bw() +
theme(panel.grid.major = element_line(colour = "transparent")) +
scale_fill_continuous(name = "Anteil ü80 [%]") +
theme(legend.position = c(0.87, 0.85)) +
theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
legend.key.height = unit(0.5, 'cm'), #change legend key height
legend.key.width = unit(0.5, 'cm'), #change legend key width
legend.title = element_text(size=9), #change legend title font size
legend.text = element_text(size=7)) +
coord_sf(datum = NA) +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "bl", which_north = "true",
height = unit(1, "cm"), width = unit(1, "cm"),
pad_x = unit(0.1, "in"), pad_y = unit(4.2, "in"),
style = north_arrow_fancy_orienteering)